{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Creating Models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One of the main goals of BIDMat/BIDMach is to make model creation, customization and experimentation much easier. \n",
    "\n",
    "To that end is has reusable classes that cover the elements of Learning:\n",
    "\n",
    "* Model: The core class for a learning algorithm, and often the only one you need to implement.\n",
    "* DataSource: A source of data, like an in-memory matrix, a set of files (possibly on HDFS) or a data iterator (for Spark).\n",
    "* DataSink: A target for data such as predictions, like an in-memory matrix, a set of files, or an iterator. \n",
    "* Updaters: Update a model using minibatch update from a Model class. Includes SGD, ADAGRAD, Monte-Carlo updates, and multiplicative updates. \n",
    "* Mixins: Secondary Loss functions that are added to the global gradient. Includes L1 and L2 regularizers, cluster quality metrics, factor model metrics. \n",
    "* Learner: Combines the classes above and provides high-level control over the learning process: iterations, stop/start/resume\n",
    "\n",
    "When creating a new model, its often only necessary to creat a new model class. We recently needed a scalable SVD (Singular Value Decomposition) for some student projects. Lets walk through creating this from scratch. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Scalable SVD\n",
    "\n",
    "This model works like the previous example of in-memory SVD for a matrix M. The singular values of M are the eigenvalues of M M^T so we do subspace iteration: \n",
    "\n",
    "$$P = M M^T Q$$\n",
    "$$(Q,R) = QR(P)$$\n",
    "\n",
    "But now we want to deal with an M which is too big to fit in memory. In the minibatch context, we can write M as a horizontal concatenation of mini-batches (this assumes data samples are columns of M and features are rows):\n",
    "\n",
    "$$M = M_1 M_2 \\cdots M_n$$\n",
    "\n",
    "and then $$P = \\sum_{i=1}^n M_i M_i^T Q$$\n",
    "\n",
    "so we can compute $P$ by operating only on the minibatches $M_i$. We need to be able to fit $P$ and $Q$ in memory, their size is only $k~ |F|$ where $k$ is the SVD dimension and $F$ is the feature set. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model Class"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We start by defining a new model class which extends BIDMach's Model class. It will always take an Options instance as an argument:\n",
    "\n",
    "<b>\n",
    "<code style=\"color:blue\">\n",
    "class SVD(opts:SVD.Opts = new SVD.Options) extends Model(opts)\n",
    "</code>\n",
    "</b>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The options are defined in the \"Object\" associated with the class. In Scala \"Object\" defines a singleton which holds all of the static methods of the class. It looks like this:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b><code style=\"color:blue\">\n",
    "object SVD  {\n",
    "  trait Opts extends Model.Opts {\n",
    "    var deliciousness = 3\n",
    "  }\n",
    "  \n",
    "  class Options extends Opts {}\n",
    "  ...\n",
    "</code></b>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Truthfully, an SVD model doesnt need a \"deliciousness\" option, in fact it doesnt need any Options at all - or rather what it needs is inherited from its parent. But its there to show how options are created. The Opts are defined as a trait rather than a class so they can be mixed in with the Options of other learning classes. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Local Variables and Initialization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are three variables we need to keep track of:\n",
    "\n",
    "<b><code style=\"color:blue\">\n",
    "  var Q:Mat = null;                                        // (Left) Singular vectors\n",
    "  var SV:Mat = null;                                       // Singular values\n",
    "  var P:Mat = null;                                        // P (accumulator)\n",
    "</code></b>\n",
    "\n",
    "and an initialization routine sets these to appropriate values."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Minibatch Update"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each update should update the stable model: Here its $P$:\n",
    "\n",
    "<b><code style=\"color:blue\">\n",
    "  def dobatch(mats:Array[Mat], ipass:Int, pos:Long):Unit = {\n",
    "    val M = mats(0);\n",
    "    P ~ P + (Q.t &ast; M &ast;&#94; M).t                               // Compute P = M &ast; M&#94;t &ast; Q efficiently\n",
    "  }\n",
    "</code></b>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Score Batches"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The score method should return a floating point vector of scores for this minibatch.\n",
    "\n",
    "<b><code style=\"color:blue\">\n",
    "  def evalbatch(mat:Array[Mat], ipass:Int, pos:Long):FMat = {\n",
    "    SV ~ P ∙ Q;                                            // Estimate the singular values\n",
    "    val diff = (P / SV) - Q;                               // residual\n",
    "    row(-(math.sqrt(norm(diff) / diff.length)));           // return the norm of the residual\n",
    "  }\n",
    "</code></b>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Update the Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At the end of a pass over the data, we update $Q$. Not all models need this step, and minibatch algorithms typically dont have it. \n",
    "\n",
    "<b><code style=\"color:blue\">\n",
    "  override def updatePass(ipass:Int) = {   \n",
    "    QRdecompt(P, Q, null);                                 // Basic subspace iteration\n",
    "    P.clear;                                               // Clear P for the next pass\n",
    "  }\n",
    "</code></b>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Convenience Functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We're done defining the SVD model. We can run it now, but to make that easier we'll define a couple of convenience functions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### An in-memory Learner"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b><code style=\"color:blue\">\n",
    "  class MatOptions extends Learner.Options with SVD.Opts with MatSource.Opts with Batch.Opts\n",
    "  \n",
    "  def learner(mat:Mat):(Learner, MatOptions) = { \n",
    "    val opts = new MatOptions;\n",
    "    opts.batchSize = math.min(100000, mat.ncols/30 + 1)\n",
    "  \tval nn = new Learner(\n",
    "  \t    new MatSource(Array(mat), opts), \n",
    "  \t\t\tnew SVD(opts), \n",
    "  \t\t\tnull,\n",
    "  \t\t\tnew Batch(opts), \n",
    "  \t\t\tnull,\n",
    "  \t\t\topts)\n",
    "    (nn, opts)\n",
    "  }\n",
    "</code></b>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A File-based Learner"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b><code style=\"color:blue\">\n",
    "  class FileOptions extends Learner.Options with SVD.Opts with FileSource.Opts with Batch.Opts\n",
    "  \n",
    "  def learner(fnames:String):(Learner, FileOptions) = { \n",
    "    val opts = new FileOptions;\n",
    "    opts.batchSize = 10000;\n",
    "    opts.fnames = List(FileSource.simpleEnum(fnames, 1, 0));\n",
    "    implicit val threads = threadPool(4);\n",
    "  \tval nn = new Learner(\n",
    "  \t    new FileSource(opts), \n",
    "  \t\t\tnew SVD(opts), \n",
    "  \t\t\tnull,\n",
    "  \t\t\tnew Batch(opts), \n",
    "  \t\t\tnull,\n",
    "  \t\t\topts)\n",
    "    (nn, opts)\n",
    "  }\n",
    "</code></b>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A Predictor"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A predictor is a Learner which runs an existing model over a DataSource and outputs to a DataSink. For SVD, the predictor outputs the right singular vectors, which may be too large to fit in memory. Here's a memory-to-memory predictor:\n",
    "\n",
    "<b><code style=\"color:blue\">\n",
    " class PredOptions extends Learner.Options with SVD.Opts with MatSource.Opts with MatSink.Opts;\n",
    "  \n",
    "  // This function constructs a predictor from an existing model \n",
    "  def predictor(model:Model, mat1:Mat):(Learner, PredOptions) = {\n",
    "    val nopts = new PredOptions;\n",
    "    nopts.batchSize = math.min(10000, mat1.ncols/30 + 1)\n",
    "    nopts.dim = model.opts.dim;\n",
    "    val newmod = new SVD(nopts);\n",
    "    newmod.refresh = false\n",
    "    model.copyTo(newmod)\n",
    "    val nn = new Learner(\n",
    "        new MatSource(Array(mat1), nopts), \n",
    "        newmod, \n",
    "        null,\n",
    "        null,\n",
    "        new MatSink(nopts),\n",
    "        nopts)\n",
    "    (nn, nopts)\n",
    "  }\n",
    "</code></b>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Testing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now lets try it out! First we initialize BIDMach as before."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 CUDA device found, CUDA version 7.0\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(0.95650464,11554000896,12079398912)"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import BIDMat.{CMat,CSMat,DMat,Dict,IDict,FMat,FND,GDMat,GMat,GIMat,GLMat,GSDMat,GSMat,\n",
    "               HMat,IMat,Image,LMat,Mat,SMat,SBMat,SDMat}\n",
    "import BIDMat.MatFunctions._\n",
    "import BIDMat.SciFunctions._\n",
    "import BIDMat.Solvers._\n",
    "import BIDMat.JPlotting._\n",
    "import BIDMach.Learner\n",
    "import BIDMach.models.{FM,GLM,KMeans,KMeansw,ICA,LDA,LDAgibbs,Model,NMF,RandomForest,SFA,SVD}\n",
    "import BIDMach.datasources.{DataSource,MatSource,FileSource,SFileSource}\n",
    "import BIDMach.mixins.{CosineSim,Perplexity,Top,L1Regularizer,L2Regularizer}\n",
    "import BIDMach.updaters.{ADAGrad,Batch,BatchNorm,IncMult,IncNorm,Telescoping}\n",
    "import BIDMach.causal.{IPTW}\n",
    "\n",
    "Mat.checkMKL\n",
    "Mat.checkCUDA\n",
    "Mat.setInline\n",
    "if (Mat.hasCUDA > 0) GPUmem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll run on the MNIST 8M (8 millon images) digit data, which is a large dataset distributed over multiple files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "val dir=\"/code/BIDMach/data/MNIST8M/parts/\"\n",
    "val (nn, opts) = SVD.learner(dir+\"data%02d.fmat.lz4\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's set some options:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "opts.nend = 10;\n",
    "opts.dim = 20;\n",
    "opts.npasses = 2;\n",
    "opts.batchSize = 20000;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and release the beast:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "pass= 0\n",
      " 4.00%, ll=-0.05831, gf=4.133, secs=0.6, GB=0.13, MB/s=206.66, GPUmem=0.948995\n",
      "25.00%, ll=-0.05782, gf=1.799, secs=9.1, GB=0.82, MB/s=89.95, GPUmem=0.948995\n",
      "48.00%, ll=-0.05741, gf=1.729, secs=17.4, GB=1.51, MB/s=86.45, GPUmem=0.948995\n",
      "70.00%, ll=-0.05785, gf=1.604, secs=27.4, GB=2.20, MB/s=80.19, GPUmem=0.948995\n",
      "91.00%, ll=-0.05776, gf=1.558, secs=37.0, GB=2.89, MB/s=77.89, GPUmem=0.948995\n",
      "100.00%, ll=-0.05784, gf=1.689, secs=37.1, GB=3.14, MB/s=84.43, GPUmem=0.948995\n",
      "pass= 1\n",
      " 4.00%, ll=-0.05776, gf=1.685, secs=38.7, GB=3.26, MB/s=84.27, GPUmem=0.948995\n",
      "25.00%, ll=-0.05781, gf=2.019, secs=39.1, GB=3.95, MB/s=100.97, GPUmem=0.948995\n",
      "48.00%, ll=-0.05746, gf=2.357, secs=39.4, GB=4.64, MB/s=117.86, GPUmem=0.948995\n",
      "70.00%, ll=-0.05781, gf=2.691, secs=39.6, GB=5.33, MB/s=134.55, GPUmem=0.948995\n",
      "91.00%, ll=-0.05788, gf=3.021, secs=39.9, GB=6.02, MB/s=151.03, GPUmem=0.948995\n",
      "100.00%, ll=-0.05750, gf=3.139, secs=40.0, GB=6.27, MB/s=156.95, GPUmem=0.948995\n",
      "Time=40.4070 secs, gflops=3.10\n"
     ]
    }
   ],
   "source": [
    "nn.train"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The model matrices for this model hold the results. They are generic matrices, so we cast them to FMats:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "val svals = FMat(nn.modelmats(1));\n",
    "val svecs = FMat(nn.modelmats(0));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfQAAAGQCAYAAABYs5LGAAA1fUlEQVR42u2de4xV1dnGSQgxxhhD0jSmaZqmSWOapt8/TfiraQjRYGxMbY0VvCIKVikUBMUbSrVa7EUQPkQsKiqlo5SbfIqKxSqi1qIDglpFkIsiNwFRVMTh/ebZhzOec+Zc9t5nX9be+/dLJrrPHGaeedfZ+1nveteljwEAAEDm6UMIAAAAMHQAAADA0AEAAABDBwAAAAwdAAAAQweA0k3Up4/vr3rvLxK/+MUvCvu3A2DoABh65tm5c6dNnTq10J0ZAAwdIGeGTnwwdAAMHSADhgUYOgCGDlAAQ2/0/iNHjtjkyZPt29/+tvXr18+OP/54GzhwoPf/jX5+o5/VTFPt97q6uuz222+37373u9a3b9+e961fv94uu+wy+853vtOj54c//KHdeOON9sknn2DoABg6AIZe7/2DBw8OPGwfhaH/5Cc/6fXeuXPneubeSIOMfd++fYxmAGDoAPk19DCT4ubNm9fLZPfu3Wu7d++O3dBrvzo7O6uuJ06c6GXxM2bM6PU6hg6AoQNg6BWv1WbnK1asCGXOYf7N/Pnzq4bQL7jggqrvHz582Htdpl75+ve+9z0MHQBDB8DQK19Tfbq2rp2UodfSv39/X3+L6uoYOgCGDpBbQw/z/rAGHIehN6udt2vIGDoAhg6Qa0OvzdDTNPQTTzyx6vuafZ9WrAAAQwfIlKEPGDCg6rW3337b18+vXdL2wgsv2Jo1a9oy9LPPPrvq+w899FCv9xw6dMhbUoehA2DoABh6xfu1/rzyNa3/Vh19+fLlTX/+D37wg8iWupV55ZVXqr7/zW9+0zo6OrxMXZpWrlzpLVvza8jsogeAoQMUxtC1plsbuLQyvcqNX8Qtt9zS6z2afd6OoQstozvuuOMiMWEMHQBDByiMoYv333/fLrroIjvhhBO8ofSf/vSn9vzzz1e9V7X2WsaNG+e9LgMePny41zlo19DFli1b7Oqrr/bKAfr56kzod/zoRz+ykSNHepk8hg6AoQOAD954440q0wtTtwYADB0AEkR7uM+cOdPbHU689dZbXpZeaegaCgcAwNABXL4ZWwxJDxs2jCABAIYO4DrnnXee/fjHP/bWgatWra9vfetbNmTIEFu2bBkBAgAMHQAAAEMHAAAADB0AAAAwdAAAAMDQAQAAMHQAAADA0AEAAABDBwAAAAwdAAAAQwcAAAAMHYLT2dlZdT1ixAi7+eabnfv67W9/66QutBEztBGzvGkbO3Yshp5F/vWvf1VdX3zxxU7qfPbZZ52NIdqIGe1JzPKkTaaOoefA0NVjdJH33nvP2RiijZjRnsQsT9ow9JwYetCGBACAfIGhk6HTy0YbMaM90UaGDq4YOjV0tBEztBGzYmvD0MnQ6WWjjZjRnmgjQwdXDJ0aOgBAscHQydDpZaONmNGeaCNDB1cMnRo62ogZ2ohZsbVh6GTo9LLRRsxoT7SRoYMrhk4NHQCg2GDoZOj0stFGzGhPtJGhgyuGTg0dbcQMbcSs2NowdDJ0etloI2a0J9rI0MEVQ6eGDgBQbDB0MnR62WgjZrQn2sjQwRVDp4aONmKGNmJWbG0YOhk6vWy0ETPaE21k6OCKoVNDBwAoNhg6GTq9bLQRM9oTbWTo4IqhU0NHGzFDGzErtjYMnQydXjbaiBntiTYydHDF0KmhAwAUGwydDJ1eNtqIGe2JNjJ0cMXQqaGjjZihjZgVWxuGToZOLxttxIz2RBsZOrhi6NTQAQCKDYZOhk4vG23EjPZEGxk6uGLo1NDRRszQRsyKrQ1DJ0Onl402YkZ7oo0MHVwxdGroAADFBkMnQ6eXjTZiRnuijQwdXDF0auhoI2ZoI2bF1oahk6HTy0YbMaM90UaGDq4YOjV0AIBig6GTodPLRhsxoz3RRoYOrhg6NXS0ETO0EbNia8PQydDpZaONmNGeaCNDB1cMnRo6AECxwdDJ0Ollo42Y0Z5oI0MHVwydGjraiBnaiFmxtWHoZOj0stFGzGhPtJGhgyuGTg0dAKDYYOhk6PSy0UbMaE+0kaGnx6ZNm+zkk0/uuV66dKmddNJJ1q9fPzvllFO86zy87tfQqaGjjZihjZgVW1smDX3KlCk2YMAA69PnazmDBg2yxYsXe/8/b948O/PMM3PxOhk62tCFNmKGttwa+pw5c6yrq6vK0E844QTvNaH/nnjiibl43a+hU0MHACg2ma6hVxp63759q76noes8vF6LhuLnzp1rs2fPrnr9vPPOswceeMBeeukl71r/deH68ccfd0pP5fWsWbOci1f5WnGjPWlP2pP2DHKdG0M//vjj62a4WX/db4ZODR1txAxtxKzY2nJj6IMHD+6ZSLZs2TKvJp2H1/0aOjV0tBEztBGzYmvLjaEvXLjQq0Prtf79+3ummIfX/Ro6NXQAgGLDOvSMQoaONnShjZihDUPPoaFTQ0cbMUMbMSu2NgydDJ1eNtqIGe2JNjJ0cMXQqaEDABQbDJ0MnV422ogZ7Yk2MnRIms7OTs/MOzo6ql4/99xzvbpO+YOm/7pwXd7O1hU9ldfaiMG1eJWvFTfak/akPWnPINcYOhk6vWy0ETPaE21k6OCKoVNDBwAoNhg6GTq9bLQRM9oTbWTo4Iqhsw4dbcQMbcSs2NowdDJ0etloI2a0J9rI0MEVQ6eGDgBQbDB0MnR62WgjZrQn2sjQwRVDp4aONmKGNmJWbG0YOhk6vWy0ETPaE21k6JA0jXaKk6G7uHMR11xzzTXX7BQHZOhoQxfaiBnayNCLZejU0NFGzNBGzIqtDUMnQ6eXjTZiRnuijQwdXDF01qEDABQbDJ0MnV422ogZ7Yk2MnRwxdCpoaONmKGNmBVbG4ZOhk4vG23EjPZEGxk6uGLo1NABAIoNhk6GTi8bbcSM9kQbGTq4YujU0NFGzNBGzIqtDUMnQ6eXjTZiRnuijQwdkoa93LnmmmuuuWYvdzJ0erJoI2ZoI2Zk6Bi6y4ZODR1txAxtxKzY2jB0MnR62WgjZrQn2sjQwRVDZx06AECxwdDJ0Ollo42Y0Z5oI0MHVwydGjraiBnaiFmxtWHoZOj0stEWTNfRo8SMDB1tGDrEZejU0CERLr/cbPFi4gCAoQMZOtqyrK3rwgvNpk8nZmToaMPQIS5Dp4aOttg5eNC+Ov98s5EjnRx2pz2JWdG1Yehk6PSy0eaP11+3zydONBs71uzdd4kZGTraMHSIw9CpoUPsLFtmdv/9ZvPmmdWcJQAAGDoEpNHhLJdeeqmThwW89NJLTumpvF68eLGzhy0obq7Fa+eNN9qmOXPM/vtfOzxmDO2Z8fbk/sxfe2LoOcnQqaGjLXYmTLA1CxaU6ueXXWa2ezcxy7g2YpYvbRh6TgydGjraYuXIEbPzz7ctGzeWru++2+yJJ4hZxrURs3xpw9BzYujU0CFWtmwxu+qqr69fecXslluICwCGDmToaMuUNn3e7rrra11ffGF20UVmhw4RMzJ0tGHoEKWhU0NHW6zMnWu2dGm1rilTzF54gZhlWBsxy5c2DJ0MnV422lrzu9+ZrV1breuZZ8ymTSNmZOhow9AhSkOnhg6xMny42YED1a/t3282bJjZV18RHwAMHcjQ0ea8tr17zUaMqK/r+uu9HeSIGRk62jB0iMjQqaGjLTZefdXs97+vr2vRotLuccQsk9qIWb60Yehk6PSy0dYcmba2e62na9s2syuvJGZk6GjD0CEqQ6eGDrFx551mq1Y1/v5vfmO2dStxAsDQgQwdbU5rGzPGbPv2xrq0pO0f/yBmZOhow9AhCkOnho62WNAGMhdcYNbV1VjXhg1m115LzDKojZjlSxuGToZOLxttjXnnHTOdgd5Ml5atXXKJ2b59xIwMHW0YOrRr6NTQIRaeftps1qzW75s+3WzFCuIFgKGDXzgPPbprzltufX1QE+KeeKJle+5eutT2TphAe3IeOu1pnIcObWbo1NDRFgs33mj25putdX32WemwFtXcac/MaCNm+dKGoefE0Kmhoy1yjh7tdaJaU1233mr273/TnhnSRszypQ1Dz4mhU0OHyPnww2CbxixfbjZzJnEDwNCBDB1tTml7+WWzP/7Rv649ezSZo5TZ055k6GjD0CGcoVNDR1vkaOLlI48E03X11WZvvUV7ZkQbMcuXNgydDJ1eNtrqM2VKr5p4S13qADz8MO1Jho42DB3CGnpkNXTNUn7mGQIMZldcYbZzZ7B/s2mTepfEDgBDh9QzdGVk555r9umn9LKLrE3trxnuNfVwX7ouv9xsxw7akwwdbRg6hDH0yGroqpuOH2+2ZEkkP446WEa1aX92rUEPo+uvfzV77DHaMwPaiFm+tGHoZOjVXHed2cqVpSxLe3TTyy6mtscfN5szJ5yuzk6zSZNoTzJ0tGHoEMbQI6mhf/JJaZhVRq5NQmp+BxSIu+8OP5fiyy81ZGR28CBxBMDQIZUM/cUXSzObxdq1ZhMm0MsuqrZrrjHbuDG8rj//OfEOYaQx03GwEZWdyNDRhqFDIEOPpIZ+zz3eQRw9qJa+bl1bP5I6WAa1aYRGZ6AfPhxe13PPlUw9qzEbPrx0glzOP2vcA/nShqGToX+Nlil98EHlLzH7/e/pZRdN27ZtZmPHtqdLw+3qZGr4PWsx+/jjUulpzBgydLSRoUPyht52DV1GLkOv5MgRs5EjzbZuJeBF4vnnzaZObf/n3HST2WuvZe/vX7OmNIdkxAiz3bv5PEBmwNDJ0EtoqF1D7rUsXmz2v/9LL7sV06ZFuuVpqnF76CGzRYva16Wla/fem732nD/f7NFHzWbMMHv6aTJ0tJGhQ7KG3nYNXZPhNCmuFm0wcsklZvv2hfqxhaiDaVh56NDSUi/XtIVB2WmDzDqQLm0uoxGehA5riSxmkyeXJoVqpOJPf8p+exb9/iyQNgydDL00CUo1Qy1bq8f995v97W/0shuhU8kU/1mz8pEBXHZZww5cYF2qxb/7bnbas6urdC+oI3vggNmwYbnej4EsmAwdUqSzs9Mz8w7t6FZj6Oo1lj9o+q/v6zfesC+uuqrx93ftsq8uvtieX7Ei3M/P+/Wdd9qBe++1Q93ZaNb/nm1a1TB8eGQ/74DOR+/+rGbl7/9Ao1TdnZDy9WF1SP77Xz7vXGfiGkMnQy/VDGs6CL3oNq2qJW30skt89llpNrdGNyLcTCW1uGmXt1tuiU5XtxlGsZ9BYjFTzVyb6pTRyFSre4MMHW1k6BClobdVQ9cmGm++2fw977xjNmpUaUgyALmvg6nOWt6MR0v8NEPaFW1h0GYqmhQXlS7VzzWEn8Bs8UhipgmglTvkaU97bYec1fYs+v1ZMG0YetEz9PJ6YT91Qu3PXW/iXJF72bffbvbCC6X//8c/zObNy3YGoNn62hQmSl3KeEOM7qQSM6091zr8Mlq62Wx+CRk62jB0iNrQQ69DX73a7I47/L33lVfMrr+e4Nd2hnSGvFi/vu4JZZlCNeMtW6L9mfrcNBnGd649a2flN1oBAoChg1MZurKnJ5/091496JTBBFhvnetetuqtymjLyNi1ZWoEu6OlEjfpPv/8pqM1oXQpLspyDx1yuz1ffbV+x2P58uq6Ohk62sjQIU5DD11D//WvS+uFg5jYH//o++25roPp5vnPf6pfmzixNBEsbW1h0PKyq6+OR5ey3HJpwtX21OS3ehPgdH/oPslaexb9/iygNgy9yBn6+++bXXllsH+jAzs0yclnJyC3veyPPiptuKMaayUPPGC2dGk2MwBNBmuxK2BoXfrZlaMZLmZNys4bbVWrCaHbt5Oho40MHeI39FA1dO1sVm+711Y88kiiW3o6yf/9X/1hWNVa/c5JcI377iv9XXGwf39km7TEgspJzSa/6fMeV2wAMHQMve0MXTO0X3op+L8r76DlY811bnvZWur3+uu9X9cOazp6M4sZgA5T0TKtuHRpQqUmDrrYnjqAqNk99O9/t3XyIBk62jB08G3ogWvo5eU42uIyDMrsFyxo+bZc1sE+/LB0ElejNfkqY1QeQ5uktnYob5ATly4d+KJthF1rT7FihZl2tWuENhDS/RJywiM1dLRh6BBfhq5MrJ0laKq/q5be4gGXy1621ps3M6bp081WrsxWBrBrl6+JX23p0vruoHM2koqZzLzVyWrah0Fb45Kho40MHeI09MA19Ci2tNTM5cpdtYrCuHFmb7/d+Pu124dmAa0VL+94Fye/+U1peNs1/Ky/V0euyS56ABg6pJOhX3NN++d3K8vXg7DJ8Zi562Xrod8qy5Rhab1+ljIAnf/997/Hr2vuXLOFC93KmlR20nB6q2NeN240u+oqMnS0kaFDvIYeqIb+8cfRzTjWuusm+5fnrg6mkY1WR8nKGBRfxTlJbe2gc7997IbWti51AjWh0JX2FFqq9rvftX6f2lUTHhscLetUe2ZcF9oyYOgDBgzofsZ9bJByhq4NPgJsDtPyZ02eXJxetrJzP1uj3nZbaRg7KxmA1ln72FugbV3qRGr9fghTjE2blmH6GJ3wmDpVT/P02lOdisWLI1v+RxZMhh6adevW2Te+8Y32jvqEuoYeqCE1Acjvdq9+HtBXXGG2aVP+g64T51Ri8IPqrQ8/nI2/SzO4L7yw9ZBzVNx1V2lWuSvcemtp21c/aLJjzBvkNEUrLM49t+1VFJBPUhlynzt3rvXr18/mRXQyFYYeMEO//PLSgyEqtOFGg4dcrnrZmtkuo/aDhpbbOKgl0bhpLoXPFQ+R6NLQ/h/+4EbWpE5MkHPstUOght0Ddn4ia091EnW+fESn15EFk6FHxve//33r27dvz5dMHsIZuu8aupYOaXg16gxPD7k9e3p9Kzd1MK0519pzvx2hNg9qSTRuOnzE585/kegqr+kun1KXRntW3g9BJzBqYpz2vU+6PfVZ0lJRja61sclNLu9PtKVn6LNmzbLjjjvOFqsWBMlm6Mqm49i2VZmDZjDntZetXeGCTubS+0OuJEg0btokqNUa7Kh1aZhbu6+lnTX52L++Fw8+GHimfiRxe/75kpFH2CEiCyZDD83LL79s/fv3t4maGQ2RGrrvhtRkre52iJy9e0uTnWI+IjM1ujuhtmxZsH+jDs6SJe7/bep4NFtXH9eoQLOd2ZJs16eeCvZv1q4tbZObNNrYpjzRUgfJNFldAsUkUUMfOHCgfdJia0mIMUPXkF2c51Jrh7Sak8Zy0cvWNrkqKah+GqwHG3qzlsTiplKCSgOff56sLpVnLr000ol4obRpk6Cg/658HylTTqo9VRrQ3JfydsPqXEYw0kYWTIYODhq6rxq6ho1vuCHOT3lp+9CKJTW5qIPpzPMwp9nphLEQE6gSjZu28B09Oh1dOns9grPjQ2srbyjTaE/+ZmjoO8CyxLbjNmdOafOfMprl3uYZ7bm5P9GGoRcyQ9eqgsqHQhxoKFC1vjz1sjWD32eNuRfa6lSm6WoGoH0E/vKXdHRpe9yIJneF0qah8yZ7KDQlYIbcVtxUK1c5q3aEKIJtdMmCydDBQUP31ZARZ0R10a5b+j15QQ/TIMuaapkxw+397rXrnd+leFGjtd8jR6b3t6tzO39+uH+rIXAZahLo81NvIygto2RiMWDo2aWzs9Mz846ag1UuvfRSbxio3HPUfyuvt61fb0cuuMDeO7YBTO33I7vevNlb1vNhd0ar65eOnbce2+9r41qrLFq9f7cemLffHvr37dUOZDNnBv73ilsS8fhMa+XXrPH9/qjb80vtvPfuu4m1Z+X1R92f013H1nOH+X1fDBtm7x+bmBZre15zje186qle31/30EP2+bHOc9j4Zf3+TOs6qfsz6DWGnpMMvWUNXcPg2q87CbSbVrcJiszXwe64w+y558L/kjDrnJOMmzJkrVBIS5c6phGdYBZIWwT77XslA587LoaOm9a7aySg3jyM8uQ8zQVIImYu3p9oI0PPo6G3rKFrrW3YOnBQ9KCRUXSbWabrYHpQqqMUYDZzQ+M4cCBabVFQPqQnTV3q8Gjr4AhmuwfSpnkN7Q6Zr17texVD6Lip09Bs6aN+v+ZBJBEz1+5PtGHoeTX0lg0pg925MzmBixZl70zwWv75z0ATxhqi0YqIN1GJBK16CDspLEq0dEz75CeJRpG0zLIdNK9CHb6IDkqp26FsNYqgTnq7fwfkBgy9CBl6khN4ymi/gUsusW0yjaz2sjVjP4pNeNS5CTisnEgG8NhjZg88kL6uBQvq7jIYa9ak3fGiOKBIm/K88UY8cVN9v9VBMCqXtLGenyyYDB0cNPSmNXQtsfnrX5MXed99tiWGQziiomkdTGvIlR2F3Iu9Cj3wfR5+4ktbVCizC/h7YtGl4W9tmtLmsHsgbdqPXRM420VzAHzMlA8VN41cvPlm6/eNHx96hIM6db60YehFyNADboIRGTt3Wpd2IZM5Zq2Xrewo6B7fjTh8OPBBLYlkADKCgKYWmy6dIBZy3/vA2srHxUYxVC7N11wTfdxk5DJ0P2jpYc2qF7JgMnQMPcOG3rAhQ2xTGSnXXVea5Z1k/T4KtKNeZ2d0P08Zup9sKym0ne3550czAhEFOuzkvvuS+V0qA0W1F7s6Be3Olq+Hhtr9HpHqs1MB+QdDz3uGrodXG+dyR9KT1cQdHfu4bl02etm7d5f0RjnZSSd0qZbuSgagzFwZuiuZiY6l1cTNNobdfWvTRjraNTEqtBx01aro4lZefeD3zAVtXaud5PbtIwsmQ8fQ82DoDWvoOtZUk45SoqfWpOxUD2wd3+oIDetg2kwm6jkHmlwXYD5B7DU6/fwQs6Nj1aUsc8OG6NuzFq06iLIEpQ5ri/JMoLhpmZpOgQua0WtVhkvtGdf9iTYMvbAZumqTSS8JatST1Qlbemjr4efAUG/DXrZiFvXwuNahK4vymYHGngFoVrlmubuUmcjI2uhI+dKm+KsdAu4L0JRdu8xGjIgmbtKnFSnaUCYI2jjqz38mCyZDx9DzYOh1G1IPLQ3dhTlNKi40QUzZhGrrIYYIY2f79tIpVhEe69mDTjXTz3cBrT93rATiGaNKHXF+XnVK2ahR0f9czRPZsqX9n6N5GxMnBv935TXxmhsB8aLnVs0x0Rg6xJ+hh+y1J9KT1bC2liq5MnpQJsKtSHuhkYkVK9zIAEJO5Ipdl9Z1h9y7wJc2DZW2WtsdBk3oa/KQ9x03HcISYujcQxMvA8aOLDgEjz5qR4cMiXaUB0PH0CupW0PXSV8+DSQumtaadNqWMrKU6lF1tSmLjmJ9cj10apbPpXCx1ui0GUnIU85irx2qDKBNX+LSpiNP/c4eD4IOadFGRO1oU7uoHKAT/sKgyX6afOlSe8b17EgLjdyNGWO7tAom6DwHDB3aytBV19OMbZd72RoC1XCldixLuDTQS5vqliEOUvGNhtt97tgXa3Yi87ntNjezJs2zGD481AoDX9o0PyJofdoPn39eWh7awIx9aXvkkfaW7qkj2upMBzL09u+da6+1rTqGWp3iuDr/GHqxDb1XQ6qeF6c5RYmW52jmsTIcbRmbFspu9FCNs3evDCztjXaUyWkzEldR9rN2bfQ/N8oNZeqheQmvvRbu30qTSlDaprmdz1fSZzYUDT2jVMoUGnGLaj8DDB1Db5qha+hyzpzs9LL1MNLaYGWw7TzUwmrT79cDVduQxolOx/KxP3ys2YkOnAl5QlciWZOWNoY42KeltvXrzSZNik+35oU02Bu/pTYd3hOFOWgYOEBJgQw9AHouqcPU3fnytOmZofPoj50pj6FDZIbeq4Z+661m//lP6joD18G0QYcOm0jgdLIqbdpvXTdn3Oih7+Mgkljrhxq5CdlxSaSu+dFHoYbdW2rTbnTalyEuNPw6dmw4bdqeucXmNL5QZ1E/y6X2TOrZETea26HPUKU27ZugVROO7LiIoecxQ9fSsDS3e223l71pU+mM7EcfjWf5WD1tmizV7NzpqND6di3ZSys7Ua1X+8qHnK+QWNakTFqTJqPUpo19ojg9rxH6rGqSp+YBBNGmIXL9uyhMQfd8k1o+GXpItCxQK0P031ptWkkUYBdIDB1aGnpVQ6r+6FhtJzBaEqIta7WtpkwoTpQJKiOs8yCOHD20ZajqdKXB22+Xloa5zvLl0R2OU0ZtHPfeB9p9T7XVIGjUIMqtaFXn1eQtiA5l5o1WX6hDps+WA4dQYeh5zNC1jloTn7Ley9YmGbqJdNRlDBN9erRpIlOS+91r0leLM7Rjy060TWnIZWGJZk0hjq9tqm3HDrMrr4xf93PPleYo+NWmv0/ZeZSfbx2XrBEnMvToOvyqnVfM7emlTR2yEPM+MHSoa+hVNXQdurFxoxM6I6mDPflk6aEXcsORltqUVel3JIU6XC2G6GKrH+pBr+w3zfb0i2aNB5gH0lSb7pepU5PriNSUNBpq04zpADVvX2gpqHY7dK09Xf6sNUNtVLPHQC9tKnU4sIwNQ89bhq4hxQB7hmeml62MVg9KPfxUh45gzbqnTUPf+rlJ7vqkg0G0TC+N7ES7ibVx7niiWZNGE+66Kxpt2iP+8ceT0a217ipt+NGmuQJxTADVapGtW8nQo0AlqpoSRl1t2uEv5VInhp4TQ+9pSL1+5535/KM1o1QPZh3wopqVdsJ78cX2Jv9pyUnUGVIrysdjJt3p0u/TOmwHJksGilMUk8W0giGpUSsNv2pCZys0hKulknGsi7///tKKCmgPbSCjVSF+7lW9R88mPZMwdIgkQw8zKSeLvWxtk6kMTpmuZvVqmZ6GkgPsjOdp0wzVlSuTD4weEk3W28cSN9WR2zyYJPGsSUOdPmemN9SmGd+aiJjUwSVa7655Eq20aVc4P8YfBk2M9bHmngy9BUqO6pSoGmrT6KHmaqS0jA1Dz4mhezX0Jstm0iKROpge2Bq21KQU/f2aQ/D3v5cOfmnSs16lDoHipp3qkkZa9fuTjJsyBx3+4Xp7VqLOqc8Rp4baNLKT5KRHdRzUyfz008ba9JlVaUxr7uNAhlKjwYn2dPmzVkt5G+I6q2yaaktxGRuGnqcMXb3GgHs5566XLQNX/XL+/NLseO1nr92zNLmqZm3uHq0EaNPgQqN6m0ZTkoybOjltZoSJt2f5WFAf66obatPDNa4T9BqhkaOKHcR6aVNHJe7PnnYlbLEjIBl6E7ScsMHnpqk2HQOcxBJJDD2/hu41pDZGUe0Mqm8uTYbS0K0yFj3kdAKdbjbVztOqd2mntjjO5W71gNeEvKzRbjvp7056e05tv9pseaDOPI9jv/pKNALUpNMITVAHUqYc9nArnZUwcyaGDm1k6A5uKOFUBqCh9dWrSzOnhw0rnWmc1gYvGklo0ouPJW5ayqQOTtbaU3McNIwZVpu2Eo5raLsRWjqm3Q7radNpb5qFHvekSM0z0d/e5PeQoTdAy1ibfOZaaisvY4vjZD8MPf+GfolmLysDjXtXtYA4W6M7csTWNzhII9GMuUHmGHncdIpd7X7/WWlP1YGlvcVnu642bdjic0125Ghy1LE986u0af5EEtsMC80n0VySrN2faWpTB0gJkia4taNNHdEk525g6Pkx9AmaxRuwMQuXobumTQ/1Bge1RK5NE8MiWCObWsy0D3uYerA2BUlrGaeG3I+tfe/Rps6JluJpSV4SaOi3o4P7MwjaOVLLz9rVVl7GplFBDB2CGPrNZ53lzAEB4BOtcU1qX3UdSaplUllFW6qGmUSmI4T1t6eBltvVbiCk2vq0aclp0CZCLcwJatASWH3eoiDhZWwYek4MfdKAAYnXa8jQ26R8UEudGdyRa9NBJxHsT5BazFST1LB7k01x6mqTmTUZco4VZeMqg3W3c4+2ceOaDuVGjnZU1PK4JOdqZPn+3L69VPtusWdBIG0aIUrobA0MPS+G/j//48x2r5VQo2uBamwaDo9bm3ZKi6DDl2rMlKE3yZzqrvVWhynNs6rVvq+/XtImI5ehJ41GBLRMkvuzNbNn+zLfQNo0Uz6hZWwYek4M/fpBg+hlZ1Gb1rrWeYBEqk1bi0Z0ZGuqMVu1qlRL96tNBqq969NE6/6729jTJmPVkHvSaB5Bgxnb3J8VlM889zG/IbA27YsR9XHAGHr26ezs9My8o2aiy/gLL/R6jeUPmv7LtfvXuzRp6rbb4v19W7bYl6NGZT5ezz/9tHVpNcenn/p6/6Y777SPj63DTkv/Do0oTJhg2zZssCPdnaqtxw7GSVRPt1Epbs91Z+ncf00+L92dnk+OzdOI+udvfftt++Lii22HOqUx/j0Yek4y9N86tkMcGUDArKCmXBKpNplKRBOxUo9Zk733e2nTwznFgzI8jtWwD2jGu3YsTAuNVNQ5fpj7s2IUSwfldHd+Y9OmYfqYl7Fh6Dkx9IsjWGMcB9TofKDOWM2DJFJt2r4yonXPqcdMBn3bbf60aV9/ba6SNn/5i32pyXFpTlpVWefBB7k/G6ElkZMnx6tNnXbtENhi+SWGjqGToWdZmzYaeeqp+LRpB8HOznzETBPd1HnVyEYzbdoRTxmXC8yda11pd7g3b657zgP35zGuu6503kPc2lRy0Q6CMe1QiaHnxNBvdnBTGfCJhpDj3HNbQ/opHBQRG1oG1GoJnrIgH9vFJsKBA6UHeZooO9RyLO2cB9XoMKcktuKt/PzGtIwNQydDJ0NPG51Trs0notamB5QO6zn//FKNMC8x04YtGnVopk1/92OP8VmrRDX8mln23J/HDDbg6oO2tJWXscVwvgCGnhNDp4aecW01B4i0rU1D03fcUaoLtjgTO3Mx03BlneVFVdq0A5924uOzVt0R0sl13J9fozkWMtcmGxbFok1HGc+YgaEDGXoutdXMyG5Lmx5S2khGM6sjysydi5lOzNPxoPW0yfDT3lDGxbjJtDQ5r2JnwsLfn/PmNTxPIVZtagPN8di4EUMHaui5Y+nSaM6y10xqPSiWLct3vHSue6NZyapXa5IT9MbBI5ZTo3zmeZtHCrfxEDe74YZIa/cYOhk6GboLaGKOlrS0o00ZvobuY3xgOxMzZd8adt+/v7c2dY7SPhrX1bipo3fvvdyfQiM8f/pTetpk5CoNHdtsBkPH0Hughp5xbToMouKglsDaNGtWE+t8boyRi5ipBrl8eW9tmt0e41rfTMftgw+qzocv7P0pMx07NvRBOZFp0zwPtUedA5owdDJ0suAsa5s0qWc3L9/alKlqyZt2AqvIVgsRs1dfrTrjvUeblmdpJjGftfpoidbWrcW+P7UvQxvHykaqbcWKyCauYug5MXRq6DlAE3QWLPD/fq1v1laSU6c6NQEsMTThT0eDVi7/kZHL0KExmquxeHGxY6DZ/jXP0DyAoZOhk6G7gmrfx5YVtdS2bZvZqFFmjzyS6LG5zsVMu+zpgJuyttWrQ9dFCxO3tWtLo0FR61KWqcxXZ4q//35pf4UPPyxtZqOOllZfqPOlTY7UGdWyQ+34p3936JDZ55+Xhp7VOT1yxLbGtexQ2kaMaHnmeRafaxh6TgydGnoOtOnhpnbsNuim2l57rbRPeQp1YudiJnPSTOGyNk2G06Q4PmuNkWFq+Vq3kUaiS+asORyakKmJijrzXfXpMWPMRo8uDfFrfodqxVqBoREUfX71fs0y17/R514n6WkeyXnnmQ0danbOOaVZ+Tr+NaIas4cmBepY2xw+1zB0MnQydJfQg7BbV0Nt2tFKD8R33iFmQsPuMoU9e0ratFwt5ESnQsVtyhSvQ9iWLq2h1vwNGfLs2T11+ajYos+4Oq1/+EPpd+h3qQOnE+zC8sknpZ+lTkgOn2sYek4MnRp6TtBmME8+Wd+4lFmMH+/chC8nYqZtXpV5KsOL6eCLXKElW2HOD1CMdRyvllsp89YyuAh3ImyIhue1okGTPzVcrs1gNm0K/nN06uDMmbltVgydDJ0M3SU0lDdtWrU2PTA19KhMJeAWlYWImVYGdBvMDsVORsNnrTWqZ196qb2nU9j8oNq3tivVULnmeWiFQcxzNxrGTLV5DZlrOF8jWgsX+uvkqlOsYf8I2oIMHWI1dGroOdGmh9UVV3ytTROK9NBSRpLg5LdMxUxDsN1Gs0Udnvvu47Pml/Hj7dWOjubvUflCh5doNYFmx+vz6VLMNCw/Z06pHq+Jfhp50LB6PTRhMsCZ51l8rmHoZOhk6K7RbU7bVSvUw1TDizV7lhOzOtx7rx1Wxhbhrlu5j9vf/mb7Va6oRSULHU87YUKpM6kSkGaguxwzZd9aJaIlnEpudDaCDqOpXM6p4XptGZzj+wBDz4mhU0PPEVp2pXq5hjePbTQDLdiwwWzIkPT25c4i2vO+cnMVxe6hh0qTDHVSX1Y/eypLKYNWmUojCzo2Vp2SJM88x9CBDB1tHtOnW5eW8GibTmLmj+6s8ks9sCM+XS7Xcevqsi5lsxqK1qx3GfnDDzsz6TKSmGnNuyZMqnOs0kHO7wMMPSeGTg09R9r277cX9RAiZmiLO5mV0V1xhdk//+nc6oDIY9bOcreMtCeGToZOho42YlZQbVvXrydmZOjgmqFTQwcAKDYYOhk6vWy0ETPaE21k6OCKoVNDRxsxQxsxK7Y2DJ0MnV422ogZ7Yk2MnRwxdCpoQMAFBsMnQydXjbaiBntiTYydHDF0Kmho42YoY2YFVsbhk6GTi8bbcSM9kQbGTq4YujU0AEAig2GToZOLxttxIz2RBsZOrhi6NTQ0UbM0EbMiq0NQydDp5eNNmJGe6KNDB1cMXRq6AAAxQZDJ0Onl402YkZ7oo0MHVwxdGroaCNmaCNmxdaGoZOh08tGGzGjPdFGhg6uGDo1dACAYoOhk6HTy0YbMaM90UaGDq4YOjV0tBEztBGzYmvD0MnQ6WWjjZjRnmgjQwdXDJ0aOgBAscHQydDpZaONmNGeaCNDB1cMnRo62ogZ2ohZsbVh6GTo9LLRRsxoT7SRoYMrhk4NHQCg2GDoZOj0stFGzGhPtJGhgyuGTg0dbcQMbcSs2NowdDJ0etloI2a0J9rI0MEVQ6eGDgBQbDB0MnR62WgjZrQn2sjQwRVDp4aONmKGNmJWbG0YOhk6vWy0ETPaE21k6OCKoVNDBwAoNrkz9E2bNtnJJ5/cc7106VI76aSTrF+/fnbKKad411l6nQwdbehCGzFDW+EMfcqUKTZgwADr0+drmYMGDbLFixd7/z9v3jw788wzM/W6X0Onho42YoY2YlZsbbky9Dlz5lhXV1eVoZ9wwgnea0L/PfHEEzP1Ohk62tCFNmKGtsIZeo/ICkPv27dv1fc0pJ2l12vp7Oz0zLyjo6OXoavXWP6g6b9cc80111wX5zr3hn788cfXzXyz8joZOtrQhTZihjYy9G4GDx7cM8Fs2bJlXq06S6/7NXRq6GgjZmgjZsXWlntDX7hwoVef1mv9+/f3zDJLr5Ohow1daCNmaCusoRcB1qEDAACGnkNDJ0NHGzFDGzEjQ8fQc2Do1NDRRszQRsyKrQ1DJ0Onl402YkZ7oo0MHVwxdGroAADFBkMnQ6eXjTZiRnuijQwdXDF0auhoI2ZoI2bF1oahk6HTy0YbMaM90UaGDq4YOjV0AIBig6GTodPLRhsxoz3RRoYOrhg6NXS0ETO0EbNia8PQydDpZaONmNGeaCNDB1cMnRo6AECxwdDJ0Ollo42Y0Z5oI0MHVwydGjraiBnaiFmxtWHoZOj0stFGzGhPtJGhgyuGTg0dAKDYYOhk6PSy0UbMaE+0kaGDK4ZODR1txAxtxKzY2jB0MnR62WgjZrQn2sjQwRVDp4YOAFBsMHQydHrZaCNmtCfayNDBFUOnho42YoY2YlZsbRg6GTq9bLQRM9oTbWTo4IqhU0MHACg2GDoZOr1stBEz2hNtZOjgiqFTQ0cbMUMbMSu2NgydDJ1eNtqIGe2JNjJ0cMXQqaEDABQbDJ0MnV422ogZ7Yk2MnRwxdCpoaONmKGNmBVbG4ZOhk4vG23EjPZEGxk6uGLo1NABAIoNhk6GTi8bbcSM9kQbGTq4YujU0NFGzNBGzIqtDUMnQ6eXjTZiRnuijQwdXDF0augAAMUGQydDp5eNNmJGe6KNDB1cMXRq6GgjZmgjZsXWhqGTodPLRhsxoz3RRoYOabF27dqq6wkTJniNyRdffPHFVzG/xo4di6EDAAAUDQwdAAAAQwcAAAAMHRqydOlSb6Kca18dHR1O6kIbMUMbMcubttq5VRh6Rqmd9Y4utBEztBGz4mnD0HNA0J5Z0XWhjZihjZjlURuGDgAAUDAwdAAAAAwdAAAAMHQAAADA0PPG0aNHbfr06Xbrrbfa5MmTva/9+/enrmvRokXeFoSnnXaanX766Z6+PXv2OBe/p556ygYOHOiUps2bN9vEiRNt8ODB3pcLqO0mTZrUo+nqq6+2HTt2pKppy5YtdsYZZzh5P9TT5so9UU+bK/dEI21p3xP1dLl4T2DoGWfevHn26KOP9lxrbeS0adNS1zVu3Dhbs2aN94D96quv7OGHH7Zf//rXTsVOs1Svv/56pwx9586dNmzYMFu3bp1TsRo+fLi3z4HaU18LFy5M9XQ/tVn5y7X7oZE2F+6JRtpcuCcaaUv7nmiky7V7AkPPASNHjrS9e/dW9RqHDh3qpNZTTz3VGS3bt2/3sqXPPvvMKUOX+axevdq5tlNWqYeWa+1Z23Yu3Q9+PldpxbCeNlfuidrf7co9UavL1XsCQ88wtcNP+oC5MkxbyWuvveZlKC6gIdgrr7zS9u3b5/vBmxQyHz3ANCSrdrzhhhucKaE8+OCD9sQTT9iuXbts5cqVdscddzj3kHXpfmj1uUrznqjV5tI9Ufu7XbknanW5ek9g6Bmm3o3nWk141apVNmLECGdq6OPHj/dqci7Ga9CgQVXDsnpgXHvttanr0gNr9OjR3gPrnHPO8WqHBw8edO7z79L90Oz3pn1P1Gpz6Z6o/d2u3BO1uly9JzB0MvTYuP/++53JMitvzHpfrhh6bXtqaC9tlL1t2rSp51p16qDnLpOhu3NP1OsIuXJP1DN0F+6JWl2u3hMYeoapHCYT6vH/6le/ckKbJtioZ53FUY600KSayhqw+NnPfpa6rnq1QRdr6C7dD/U+V67cE60+8y5l6K7cE7W6XL0nMPQMM3/+fG8mbxn9/9SpU1PX9fLLLzuVlWfF0NV+M2bM6Imd6qx333136rq0JEdDjOUMqZxlutZ2Lt0PtdpcuieyZOiu3BO1uly9JzD0DOPqOnSXh7VdNnRxzz33eBOANKyoB9mRI0dS16TPlD5b0qQv1Q3TrBc2+my5cD800ubCPeH396e9bK1WW5r3RCNdrt0TGDoAAECBwdABAAAwdAAAAMDQAQAAAEMHAAAADB0AAABDBwAAAAwdAAAAMHQAAADA0AHAGZ555hn7+c9/7u0aVg8dt6nv630AgKEDgMNoO1Btw7lgwYKq17Wfu16fPXs2QQLA0AEgC2gvbZn3iy++6F2vXLnSu9brAIChA0BG0CEdo0eP9s48X7JkifffUaNG2eHDhwkOAIYOAFni448/tqFDh3qZ+ZAhQzJzbC8Ahg4AUIEMHEMHwNABIMOUh9x1HrUmxzHkDoChA0AGuemmm7zM/Nlnn/WuV61a5V3rdQDA0AEgA5SXrWmZWiXK1PW6vg8AGDoAOEx5YxltIFMPbTjDxjIAGDoAAACGDgAAABg6AAAAYOgAAACAoQMAAGDoAAAAgKEDAAAAhg4AAAAYOgAAAIYOAAAAGDoAAABg6AAAAIChA2T5ZuzTJ7NfAIChA0CFoaMbADB0AAwd3QAYOgBgjBg6AIYOABg6AGDoAIChAwCGDoChh2TgwIEYOgCGDgBFMfRm78PQATB0AMDQAQBDB8DQ/ZrxokWLbMiQITZ48GC74YYb7ODBg1VG/dVXX9mMGTPs9NNP976mT5/uvVZ+T+UXhg6AoQOAX0OfPDn8Vx1Dnz9/vh06dMiOHj3qmfvUqVOrDH327Nm2YMEC7/v6WrJkid1zzz1k6AAYOgC0ZegbNoT/qmPolSjzPuOMM6q+d/bZZ/dk5OLIkSP2y1/+EkMHwNABoC1Dj5B6ZnzqqadWfa/ZezB0AAwdABw0dGXiZ511VtX3lI3XZujl92DoABg6ADhi6Js3b/b+X/Xxxx9/3GbOnFll1Kqhd3R09NTQVU+vrKGfeeaZtn37dgwdAEMHgDQNfdy4cXbaaad5M9g1m/3w4cNVhq7sXDPb9R596f+VpZd54oknembAY+gAGDoApGToWdQNABg6AIZegTJuDB0AQweAjBs6ugEwdADAGDF0AAwdADB0AMDQAQpo6Fn9AgAMHQAAACLg/wH8+/NJUO45kAAAAABJRU5ErkJggg==",
      "text/plain": [
       "BufferedImage@3d0e7bf5: type = 2 DirectColorModel: rmask=ff0000 gmask=ff00 bmask=ff amask=ff000000 IntegerInterleavedRaster: width = 500 height = 400 #Bands = 4 xOff = 0 yOff = 0 dataOffset[0] 0"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "semilogy(svals)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To see how well we did, we will compute the SVD directly by computing $M M^T$ and computing its eigendecomposition. Normally we can't do this because $MM^T$ is nfeats x nfeats, but for this dataset nfeats is only 784. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "..."
     ]
    }
   ],
   "source": [
    "tic\n",
    "val MMt = zeros(784,784);\n",
    "for (i <- 0 until opts.nend) {\n",
    "val Mi = loadFMat(dir+\"data%02d.fmat.lz4\" format i);\n",
    "MMt ~ Mi *^ Mi;\n",
    "print(\".\");\n",
    "} \n",
    "println;\n",
    "toc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we call an eigenvalue routine to compute the eigenvalues and eigenvectors of $MM^T$, which are the singular values and left singular vectors of $M$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "val (eval, evecs) = feig(MMt);\n",
    "val topvecs = evecs(?, 783 to 784-opts.dim by -1);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Eigenvectors have a sign ambiguity, and its common to see V or -V. So next we compute dot products between the two sets of vectors and flip signs if a dot product is negative:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "val dots = svecs ∙ topvecs;\n",
    "svecs ~ svecs ∘ (2*(dots>0) - 1);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets now look at the eigenvectors as small images, decreasing in strength from left to right. First the reference eigenvectors:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "val onerow = topvecs.view(28,28*opts.dim);\n",
    "val nc = onerow.ncols;\n",
    "val tworows = onerow(?,0->(nc/2)) on onerow(?,(nc/2)->nc)\n",
    "show((tworows.t*500+128) ⊗ ones(3,3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "val onerow = svecs.view(28,28*opts.dim);\n",
    "val nc = onerow.ncols;\n",
    "val tworows = onerow(?,0->(nc/2)) on onerow(?,(nc/2)->nc)\n",
    "show((tworows.t*500+128) ⊗ ones(3,3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Extracting Right Singular Vectors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So far, we obtained the singular values and left singular vectors from the model's modelmats array. The right singular vectors grow in size with the dataset, and in general wont fit in memory. But we can still compute them by running a predictor on the dataset. This predictor takes a parametrized input file name for the matrix $M$ and a parametrized output file name to hold the right singular vectors. A key option to set is <code>ofcols</code> the number of samples per output file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "val (pp, popts) = SVD.predictor(nn.model, dir+\"data%02d.fmat.lz4\", dir+\"rSingVectors%02d.fmat.lz4\")\n",
    "popts.ofcols = 100000                  // number of columns per file, here the same as the input files\n",
    "popts.nend = 10                        // Number of input files to process"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pp.predict"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Notes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "BIDMach's actual SVD model is only slightly different from the code presented here: First it computes the subspace updates on *minibatches* of data for the first few iterations to more rapidly get near a good model. Secondly it alternates subspace iteration with Rayleigh-Ritz iterations for faster convergence."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "scala",
   "version": "2.11.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
